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ABSTRACT 

We present a method to map the artificial sky brightness across large territories in 
astronomical photometric bands with a resolution of approximately 1 km. This is useful 
to quantify the situation of night sky pollution, to recognize potential astronomical 
sites and to allow future monitoring of trends. The artificial sky brightness present 
in the chosen direction at a given position on the Earth's surface is obtained by the 
integration of the contributions produced by every surface area in the surrounding. 
Each contribution is computed based on detailed models for the propagation in the 
atmosphere of the upward light flux emitted by the area. The light flux is measured 
with top of atmosphere radiometric observations made by the Defense Meteorological 
Satellite Program (DMSP) Operational Linescan System. 

We applied the described method to Europe obtaining the maps of artificial sky 
brightness in V and B bands. 

Key words: atmospheric effects - site testing - scattering - light pollution 



1 INTRODUCTION 

The night sky is a world heritage. In recent decades there 
has been a rapid increase in the brightness of the night sky in 
nearly all countries. This has degraded astronomical viewing 
conditions. The increase in night sky brightness is one of the 
most noticeable effects of light pollution, which can be de- 
fined as the alteration of natural light levels in the outdoor 
environment due to artificial light sources. The widespread 
use of artificial lighting with little regard to fixture shield- 
ing or energy conservation is the primary source of anthro- 
pogenic light contributing to the brightness of the night sky. 
The astronomical community has expressed its concern over 
the growth of the sky brightness in a number of official doc- 
uments and positions (e.g. the Resolutions of the General 
Assemblies of the International Astronomical Union (IAU) 
XVI/9 1976, XIX/B6 1985, XX/A2 1988, XXIII/A1 1997 
and the Positions of the American Astronomical Society). 
The UNESCO, the United Nations and the Commission In- 
ternationale de l'Eclairage (CIE) give consideration to this 
concern. The Commission 50 of IAU (The protection of ex- 
isting and potential astronomical sites) is working actively in 
order to preserve the astronomical sky and many studies and 
meetings have been dedicated to this topic (e.g. Crawford 
1991; Kovalevsky 1992; McNally 1994; Isobe & Hirayama 
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1998; Cinzano 2000a; see also Cinzano 1994 for a large ref- 
erence list). Laws, bills, standard rules, ordinances, regu- 
lations limit in many countries the direct wasting of light 
in the sky from lighting fixtures and, in some cases, also 
the quantity of light reflected by lighted surfaces. Lhe In- 
ternational Dark-Sky Association is active worldwide in this 
battle of culture and intelligence with the aim of building 
awareness and saving for mankind the possibility of feeling 
part of the universe. 

An effective battle against light pollution requires the 
knowledge of the condition of the night sky across large ter- 
ritories, recognition of vulnerable areas, the determination 
of the growth trends and the identification of most polluting 
cities. Lherefore a method to map the artificial sky bright- 
ness across large territories is required. This is also useful in 
order to recognize areas with low level of light pollution and 
potential astronomical sites. In the past, mappings of sky 
brightness for extended areas has been performed based on 
population density data with some simple modelling. Exam- 
ples include the works of W alker (1970, 1973) in Cali fornia 
and Arizona, Albers (Il998|) in USA, Bertiauet al. (jl973j) 



in Italy, Berry (1976) and Pike & Berry (1976) in Ontario 



(Canada). These authors used population data of cities to 
estimate their upward light emission and a variety of empir- 
ical or semi-empirical propagation laws for light pollution 
in order to compute the sky brightness produced by them. 
Recently advances in the availability and gain control of the 
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DMSP satellite now provides direct measurements of the 
upward light emission from nearly the entire surface of the 
Earth (60° South to 72° North). The first global DMSP im- 
age of the Earth at night was produced at 10 km resolution 
(Sullivan 1989, 1991) using DMSP film strips, the only data 
available at that time. Beginning in 1992, the availability 
of digital DMSP data has enabled a new set of higher spa- 
tial resolution products (Elvidge et al. 1997a, 1997b, 1997c). 
These data have been used to model population distribu- 
tion (Sutton et al. 1997) and urban sprawl impacts on food 
production (Imhoff et al. 1997a, 1997b). The data have also 
been used to document light pollution as expressed in the in- 
crease in the upward flux of light over time for Japan (Isobe 
& Hamamura 1998). 

In this paper we present a method to map the artificial 
sky brightness across large territories. In order to bypass er- 
rors arising when using population data to estimate upward 
flux, we construct the maps based on direct measuraments of 
the upward flux as observed from space using DMSP satel- 
lite night-time images and compute the downward flux to 
the Earth's surface with detailed modelling of light pollu- 
tion propagation in the atmosphere. We also present, as an 
application, detailed maps of artificial sky brightness in Eu- 
rope in V and B astronomical photometrical bands with a 
resolution of approximately 1 km. In section || we describe 
the OLS-DMSP satellite observations, their reduction and 
analysis. We also discuss the relation upward flux versus city 
population. In section ^ we describe the mapping technique 
and in section ^| we discuss the application to the maps of 
Europe. The maps are presented in section [s] together with 
our comments. Section H contains our conclusions. 



2 OBSERVATIONS AND DATA ANALYSIS 
2.1 Satellite data 

U.S. Air Force Defense Meteorological Satellite Program 
(DMSP) satellites are in low altitude (830 km) sun- 
synchronous polar orbits with an orbital period of 101 min- 
utes. With 14 orbits per day they generate a global night 
time and day time coverage of the Earth every 24 hours 
with their main purpose to monitor the distribution of 
clouds and to assess navigation conditions. The U.S. De- 
partment of Commerce, NOAA National Geophysical Data 
Center (NGDC) serves as the archive for the DMSP pro- 
gram. The digital archive was initiated in 1992. Starting in 
1994, NGDC has embarked on the development of night time 
lights processing algorithms and products from the DMSP- 
OLS (Elvidge et al. 1997a, 1997b, 1997c). 

The Operational Linescan System (OLS) is an oscillat- 
ing scan radiometer with low-light visible and thermal in- 
frared (TIR) imaging capabilities which first flew on DMSP 
satellites in 1976. At night the OLS uses a Photo Multiplier 
Tube (PMT) to intensify the visible band signal. The pur- 
pose of this intensification is to observe clouds illuminated 
by moonlight. The PMT data have a broad spectral response 
from 440 to 940 nm (485 - 765 nm FWHM) with highest sen- 
sitivity in the 500 to 650 nm region (see Fig. |l|). This cov- 
ers the range for primary emissions from the most widely 
used lamps for external lighting: Mercury Vapour (545 nm 
and 575 nm), High Pressure Sodium (from 540 nm to 630 
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Figure 1. Spectral sensitivity of OLS-PMT detector of DMSP 
satellite F12 . 



nm), Low Pressure Sodium (589 nm). The sensitivity of the 
PMT, combined with the OLS-VDGA (variable digital gain 
amplifier) and fixed gain OLS pre- and post- amplifiers al- 
lows measurement of radiances down to 10 -10 W cm -2 sr" 1 
^m" 1 (Elvidge et al. 1999). This implies that the OLS-PMT 
could detect radiation with an effective wavelength near 550 
nm down to a luminance of approximately 0.2 mcd m -2 . 
The TIR detector is sensitive to radiation from 10.0 to 13.4 
^m (10.3 - 12.9 FWHM). 

The OLS acquires swaths of data that are 3000 km wide. 
The OLS sinusoidal scan maintains a nearly constant along- 
track pixel-to-pixel ground sample distance (GSD) of 0.56 
km. Likewise, the electronic sampling of the signal from the 
individual scan lines maintains a GSD of 0.56 km. The ef- 
fective instantaneous field of view (EIFOV) is larger than 
the 0.56 km GSD and varies with scan angle. At nadir the 
EIFOV is 2.2 km and expands to 5.4 km at the edge of scan 
1500 km out from nadir. Most of the data received by NGDC 
has been "smoothed" by on-board averaging of 5 x 5 pixel 
blocks, yielding data with a GSD o f 2.8 km. Other details 



of the OLS are described by Lieske (1981). 

DMSP platforms are stabilized using four gyroscopes 
(three axis stabilization), and platform orientation is ad- 
justed using a star mapper, an Earth limb sensor, and a so- 
lar detector (Elvidge et al. 1997a) . Daily radar bevel vector 
sightings of the satellites provided by Naval Space Command 
allow Air Force orbital mechanics models to compute geode- 
tic subtrack of each orbit giving the position of satellites 
every 0.4208 seconds. This position together with OLS scan 
angle equations, an oblate ellipsoid Earth sea level model 
and digital terrain elevations from U.S. Geological Survey, 
EROS Data center, allow geolocation algorithms to estimate 
latitude and longitude for each pixel center. 

Normally the OLS is operated at high gain settings 
for cloud detection. City lights are saturated in these data 
and radiances can not be extracted. In 1996 and 1997, 
NGDC made special requests to the Air Force for collec- 
tion of OLS PMT data at reduced gain settings. This re- 
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quest was granted for the darkest nights of lunar cycles in 
March of 1996 (8 nights) and January and February (10 
nights each) of 1997. During these experimental data col- 
lections on board algorithms which adjust the visible band 
gain were disabled. The major on-board algorithm which 
effects the night time visible band data is the along scan 
gain control (ASGC), which adjust the gain dynamically in 
response to scene brightness. There is also a on-board bi- 
directional reflectance (BRDF) algorithm designed to reduce 
the brightness of the image "hot-spot" which occurs where 
the observation angles matches the illumination angle. The 
BRDF algorithm has minimal effect when lunar illumination 
is low, as was the case during the 28 nights in which the gain 
controlled data were acquired. The two on-board gain con- 
trol algorithms were turned off to simplify the retrieval of 
radiances from the special data collections. 

During the special data acquisition the OLS VDGA 
gain, which is normally operated at 60 dB, was reduced to 
avoid saturation in urban centres. On one set of nights the 
gain was operated at a setting of 24 dB. This produced data 
that avoided saturation on major on major urban centres, 
but did not permit detection of city edges and lighting in 
smaller towns. To overcome this dynamic range limitation, 
data was also acquired at gain settings of 40 and 50 dB. 

2.2 Data reduction 

These data were used to assemble a cloud-free composite im- 
age calibrated to top-of-atmosphere (TOA) radiances. The 
composition provides additional advantages in the removal 
of ephemeral lights sources, such as fire and lightning, plus 
the retrieval of lights from small towns that are near the de- 
tection limits of the sensor and processing algorithms. The 
primary processing steps include: 1) establishment of a refer- 
ence grid with finer spatial resolution than the input imagery 
using the one kilometre equal area Interrupted Homolosine 
Projection (Goode 1925; Steinwand 1993) developed for the 
NASA-USGS Global 1km AVHRR project; 2) identification 
of the cloud free section of each orbit based on OLS-TIR 
data; 3) identification of lights, removal of noise and solar 
glare, cleaning of defective scan lines and cosmic rays; 4) 
projection of the lights from cloud-free areas from each or- 
bit into the reference grid; 5) calibration to radiance units 
using prior to launch calibration of digital number for given 
input telescope illuminance and VDGA gain settings in sim- 
ulated space conditions; 6) tallying of the total number of 
light detections in each grid cell and calculation of the av- 
erage radiance value; 7) filtering images based on frequency 
of detection to remove ephemeral events. 

The final image was transformed into latitude/longitude 
projection with 30" x 30" pixel size with data in 8-bit byte 
format and power law scaling (see Elvidge et al. 1999 for 
details). The maps of Europe described below were obtained 
from an image extracted from the global image. The image 
of Europe is 4800x4800 pixels in size, starting at longitude 
10° 30' west and latitude 72° north. Orbits involved are listed 
in appendix 

Radiance range of final composite image goes from a 
minimum of 1.54 x 10~ 9 to a maximum of 3.17 x 10~ 7 W 
cm -2 sr _1 /Ltm -1 . The minimum luminance detectable for 
light with effective wavelength of 550 nm (which power is 
1.47x 10 -3 W lm -1 ) is approximately 3 mcd m -2 . Assuming 
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Figure 2. Normalized upward flux versus city population rela- 
tionship for Italy (open triangles), France (filled triangles), Ger- 
many (filled circles), Spain (filled squares), Greece (open squares). 



an average vertical extinction of ~0.3 mag in visual band, 
the minimum detectable luminance on the ground would be 
of the order of 4 kcd/km 2 . Two un-shielded fixtures of the 
"globe" kind, with clean transparent glass (fixture efficiency 
~80 per cent), equipped with a 250 W high pressure sodium 
lamp with 125 lm W _1 efficiency, and placed every square 
kilometer, could be sufficient to produce this luminance. 



2.3 Data analysis 

We first analysed the composite image in order to evalu- 
ate the emission versus population relationship. We chose a 
number of European cities of various populations^ and mea- 
sured their relative upward flux per unit solid angle summing 
the counts of all pixels pertinent to each city and multiply- 
ing for pixel size at that latitude. Fig. ^ show the measured 
upward flux of a sample of European cities normalized, for 
display purposes, to the average flux of a city of ~ 10 inhab- 
itants in the same country. The upward emission increases 
linearly with the population. One possible source of error 
is the on-board averaging of the 5x5 pixel blocks during 
the smoothing process. During smoothing it is possible for 
saturated pixels in the cores of urban centres to be aver- 
aged in with unsaturated pixels to produce an unsaturated 
smoothed pixel value. This phenomenon will be addressed 
through the inclusion of OLS data acquired at even lower 
gain settings in the updated nighttime lights map NGDC 
is preparing for the 1999-2000 time period. Another uncer- 
tainty in the city analysis is that we did not precisely match 



T Population data refer to 1991 for Italy, Spain and Greece, 1990 
for France, 1996 for Germany and have been provided by their 
national bureaux of statistic. 
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the outline of the lights to match the population data re- 
porting area. 

Even if the upward emission vs. city population rela- 
tionship depends on the local lighting habits, our results 
are in agreement with relations success fully a pplied or mea- 
sured in other countries. Elvidge et al. ( 199S| ) found a linear 
relation between composite radiance and population for a 
sample of cities of North America and obtained a smaller 
scatter when c orrelating radiance with electric consumption. 
Walker (1977) found linear proportionality between popula- 



tion and street light emission for a number of California 
cities, with a few departures above or below the mean de- 
pending on the industrial or residential character of the city. 
He also measured the sky illumination produced by three 
cities, obtaining a depende nce o n the power of 0.8 of their 
population. Bertiau et al. ( 1973 ) used successfully a linear 
model with brightness proportional to population to predict 
light pollution, and so did Walker (1970, 1973) and Garstang 
(1987, 1988, 1989a, 1989b, 1989c, 1991a, 1991b, 1991c, 1992, 
1993, 2000) who o btaine d a good fit with many observations 
including Walker (1977) population - distance data. Berry 
(1976) fitted well observations of sky brightness in city cen- 
ters in Ontario with a propagation law f or light pollution 
based on the approach of Treanor (1973) but assuming a 
dependence on the power of 0.5 of t he po pulation. However 
Garstang's linear models fit Berry's (1976) observations well, 
suggesting that the power of 0.5 found by him was produced 
by the extinction of light emitted by outskirts of large cities 
in propagating to the center and does not depend on the up- 
ward flux versus p opulat ion relationship (Garstang 1989a). 

Bertiau et al. (1973) in the early 1970's found that the 
city upward emission in Italy depended on its economic and 
commercial development, so they were forced to include in 
their model a development factor. Twentyfive years later we 
are unable to identify an affluence effect on the brightness 
of italian cities. Cities in southern Italy have the same light 
output as comparable size cities in northern Italy, even if 
the former have a per capita income that is nearly half than 
the latter. 

The proportionality between satellite data and popula- 
tion would allow us to replace the population distribution 
from census data, usually adopted in computations of night 
sky brightness, with satellite data, independently by the fact 
that we would be observing the light coming from the exter- 
nal night-time lighting or whatever else. This replacement 
constitutes an improvement because census data are based 
on city lists which, though they can be associate with the 
city geographical position, do not provide spatially explicit 
detail of the population geographical distribution. 



3 MAPPING TECHNIQUE 

Scattering from atmospheric particles and molecules spreads 
the light emitted upward by the sources. If e(x, y) is the 
upward emission per unit area in (a:, y), the total artificial 
sky brightness in a given direction of the sky in a site in 
(x', y') is: 



where f((x, y), (x' , y')) give the artificial sky brightness per 
unit of upward light emission produced by the unitary area 
in (x, y) in the site in (x' , y'). The light pollution prop- 
agation function / depends in general on the geometrical 
disposition (altitude of the site and the area, and their mu- 
tual distance), on the atmospheric distribution of molecules 
and aerosols and their optical characteristics in the chosen 
photometric band, on the shape of the emission function of 
the source and on the direction of the sky observed. In some 
works this function has been approximated with a variety of 
semi-empirical propagation law like Treanor Law (Treanor 
1973; Falchi & Cinzano 2000; Cinzano et al. 2000), Walker 
Law (Walker 1973; Joseph et al. 1991), Berry Law (Berry 
1976; Pike 1976), Garstang Law (Garstang 1991b). 

We obtained the propagation function f((x,y), (x' , y')) 
for each couple of points (x,y) and (x' , y') with detailed 
models for the light propagation in the atmosphere based 
on the modelling technique introduced and developed by 
Garstang (1986, 1987, 1988, 1989a, 1989b, 1989c, 1991a, 
1991b, 1991c, 1992, 1993, 2000) and also applied by Cin- 
zano (2000b, 2000c, 2000d). The models assume Rayleigh 
scattering by molecules and Mie scattering by aerosols and 
take into account extinction along light paths and Earth 
curvature. These models allow associating the predictions 
with well-defined parameters related to the aerosol content, 
so the atmospheric conditions at which predictions refer can 
be well known. Here we will describe only the main lines 
of the models and our specific implementation, leaving the 
reader to the cited papers for details. 

A telescope of area 7rd 2 /4 situated in the observing 
site O collects from within an infinitesimal section dV = 
(ire 2 u 2 du) of the cone of angle 2e around the line-of-sight 
at a distance u and with thickness du, a luminous flux d<I> 
given by: 



(2) 



where M s (u) is the luminous flux scattered in unit solid 
angle toward the observer from particles of aerosol and 
molecules inside unit volume of atmosphere at the distance 
u along the line of sight, and £i{u) is the extinction of the 
light along its path to the telescope. 

Calling e s the upward flux of the source and S = M s /e s 
the scattered flux per unit solid angle per unit upward flux, 
the propagation function /, expressed as total flux per unit 
area of the telescope per unit solid angle per unit total up- 
ward light emission, is found integrating eq. (Q) from the 
site to infinity: 



/= / S(u)^(u)du, 

J UQ 



with: 



6 



exp 



(N m (h)a m + N^h)^) dx 



(3) 



(4) 



b(x', y) 



e{x,y)f{(x,y),(x', y')) dx dy , 



(1) 



where iV m (/i) and N a (h) are respectively the vertical num- 
ber densities of molecules and aerosols, a m and a a are their 
scattering cross sections. The altitude h depends on the in- 
tegration variable x, on the zenith distance and azimuth of 
the line of sight, on the altitudes of site and source, and on 
their distance. 
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The luminous flux per unit solid angle per unit upward 
flux coming directly from the source and scattered toward 
the observer from unit volume along the line of sight is: 



S d {u) = (N m (h)a m f m (uj) + JV a (/iW.(w)) s ) 



(5) 



where i(ip, s) is the direct illuminance per unit flux produced 
by each source on each infinitesimal volume of atmosphere 
along the line-of-sight of an observer and / m (cj), f&(u>) are 
the normalized angular scattering functions of molecular and 
aerosol scattering respectively. The scattering angle ui, the 
emission angle ip, the distance s of the section from the 
source and the altitude h of it, depend on the altitudes of 
the site and the source, their distance, the zenith distance 
and the azimuth of the line-of-sight and the distance it along 
the line of sight, through some geometry. 

If I(ip) is the normalized light flux per unit solid angle 
emitted by the considered source at the zenith distance ip 
and s is the distance between the source and the considered 
infinitesimal volume of atmosphere, the illuminance per unit 
flux is: 



i(iP,s)=I(i>)&/s 2 



(6) 



in the range where there is no shielding by Earth curvature 
and zero elsewhere. The extinction £2 along the path is: 



£2 = exp 



(N m (h)a m + N a (h)<Ta.) dx 



A single scattering model is not sufficient to describe 
the artificial sky brightness produced by a source. In a real 
atmosphere several scatterings may occur during the travel 
of a photon from the source to the telescope. The optical 
thickness r = J kdr, where k is an attenuation coefficient, 
determines how important secondary and higher scattering 
is. If r >> 1 (thick layer) multiple scattering is dominant. 
The fraction of incident radiation which has been scattered 
once is (1 — e _T ) and the fraction which it is scattered again 



is of order (1 — e 



If (1 - 



is sufficiently small, which 



happens when r is small, secondary and higher order scatter- 
ing can be neglected. In absence of aerosol the optical thick- 
ness of the atmosphere at wavelength of 550 nm is about 0.1 
(Twomey 1977). The aerosol optical thickness can be 0.05 in 
cleaner regions of the globe, but it can grow to higher values, 
even depending on seasonal changes (Garstang 1988). Then 
single scattering is the major contributor to scattered radia- 
tion but secondary scattering is not negligible. The error in 
neglecting third and higher order scattering can be signifi- 
cant for optical thickness higher than about 0.5. Therefore 
in the computation of the total light flux S that molecules 
and aerosols in the infinitesimal volume scatter toward the 
observer we take into account light already scattered once. 
We assumed, as Garstang (1984, 1986), that the light com- 
ing to the considered infinitesimal volume along the line of 
sight after a scatter has approximately a direction near that 
of the direct light, so that the scattering angle ui in first 
approximation is always the same. In this case the total il- 
luminance S can be written: 



S = Sd Ds . 



(8) 



where Ds is a correction factor which take into account 
the illuminance due at light already scattered once from 
molecules and aeros ols wh ich can be evaluated with the ap- 
proach of Treanor ( 1973 ) as extended by Garstang (1984, 
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1986). Details on assumptions can be found in the quoted 
papers. 



4 APPLICATION 

In practice, we divided the surface of Europe into land areas 
with the same positions and dimensions as projections on 
the Earth of the pixels of the satellite image. We assumed 
each land area be source of light pollution with an upward 
emission e XlV proportional to the radiance measured in the 
corresponding pixel multiplied by the surface area. The total 
artificial sky brightness at the centre of each area, given by 
the expression (fil), becomes: 



(9) 



We obtained the propagation function f({xi,yj),(xh,Vk)) 
for each couple of points (xi,yj) and (xh,yk), the positions 
of the observing site and the polluting area, from eq. (^|) 
after inserting eq. (^) and eq. (g) multiplied for the double 
scattering factor Ds computed as below. We considered ev- 
ery land area as a point source located in its centre except 
when i = h, j = k in which case we used a four points ap- 
proximation (Abramowitz & Stegun 1964). The resolution 
of the maps, depending on results from an integration over 
a large zone, is greater than resolution of the original im- 
ages and is generally of the order of the distance between 
two pixel centres. However where sky brightness is domi- 
nated by contributions of nearest land areas, effects of the 
resolution of the original image could became relevant. 

We obtained the maps for B and V photometric astro- 
nomical bands (Johnson 1955). 



4.1 Atmospheric model 

4-1.1 Molecular atmosphere vertical distribution 

We assumed the atmosphere in hydrostatic equilibrium un- 
der the gravitational force. Neglecting the Earth curvature, 
the force per unit surface supporting a molecular layer of 
thickness dh is dp = —gpdh where g is the gravitational ac- 
celeration and p the density of the layer. Replacing dp with 
the equation of state of perfect gas, working for dry air, 
dp — dpRT/M, where R is the gas constant, M the mass 
of a mole of dry air and T the average temperature, and 
integrating, we find that in first approximation the num- 
ber density N m of the gaseous component of atmosphere 
decreases exponentially with altitude h. If n is the average 
number of molecules in a mole of dry air: 



N m (h) = ^ exp (-Mih 
M V RT 



m,o e 



(10) 



Measurements show that this is a good approximation for 
the first 10 km (se e e.g . Fig. ^). We adopted this simple 
model as Garstan g ( 1986 ), neglecting improvements d one by 
Garstang ( 1991a ) for higher altitudes. As Garstang ( 1986 ) 
we assumed an inverse scale altitude c = 0.104 km -1 and a 
molecular density at sea level Nm,o 
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Figure 3. Examples of vertical distributions of molecules 
(squares) and haze aerosols (circles) in the atmosphere from 
respec tively U.S. Standard Atmosphere (1962) and Eltermann 
(1964) Clear Standard Atmosphere. 



4.2 Angular scattering functions 

We take into account both Rayleigh scattering by molecules 
and Mie scattering by aerosols. For molecular Rayleigh scat- 
tering the angular scattering function is : 



/ m (w) = 3(l + cos 2 (u;))/167r. 



(13) 



The integrated Rayleigh scattering cross section in V band 
was assumed to be <r m = 1.136 x 10 -26 cm -2 sr _1 and in B 
band a m = 4.6 x 10~ 27 cm" 2 sr" 1 . 

The normalized angular scattering function for atmo- 
spheric haze aerosols can be measured easily with a num- 
ber of well known remote-sensing techniques like classi- 
cal searchlight probing (see e.g. Elterman 1966), modern 
bistatic lidar probing, measurements of the day-light or 
moonlight sky scattering function (see e.g. Hulburt 1951; 
Volz 1987; Krisciunas & Schaefer 1991). Nevertheless in this 
paper we are not interested in a specific function for a given 
site at a given time but in the typical average function so we 
adopted the function tabulate d by M c Clatchey et al. ( 197£ ) 
as interpolated by Garstang (1991a) and we neglected geo- 
graphical gradients: 



For < lo < 10° 
/a(o>) = 7.5 exp I 



-0.1249 lo 2 /(1 + 0.04996 lo 2 )) 



For 10° < lo < 124° 



4-1.2 Haze aerosol vertical distribution 

We are interested in average atmospheric conditions, better 
if typical and not for the particular conditions of a given 
night, so a detailed modelling of the local aerosol distribu- 
tion for a given night is beyond the s cope of this paper. As 
Garstang ( 1986j ) and Joseph et al. (1991) we assumed an 
exponential decrease of number density for the atmospheric 
haze aerosols: 



iV a (ft) = iVa.O 



(11) 



Measurements show (see e.g. Fig. |3|) that for the first 10 km 
this is a reasonable approximation. To account for presence 
of sporadic denser aeros ol layer s at various heights or at 
ground level as Garstang (1991t) is beyond the scope of this 
work. We also neglected the effects of the Ozone layer and 
the presence of volcanic dust studied by Garstang (1991a, 
1991c). 

We t ake in to account changes in aerosol content as 
Garstang (198C) introducing a parameter K which measures 
the relative importance of aerosol and molecules for scatter- 
ing light in V band: 



K = 



A'a.OO-a 



iV m ,o<Tmll.lle- cH 



(12) 



where H is the altitude over sea level of the ground level. 
Changing the parameter K we are able to compute the map 
for different aerosol contents, i.e. for different products ./V a <7 a . 



As Garstang (1986), the inverse scale altitude of aerosols 



was assumed to be a = 0.657 + 0.059ii\ Effects of ch anges 
of aerosol scale altitude has been checked in sec. 4.7.1. More 
detailed atmospheric models could be used whenever avail- 
able. 



exp (-0.07226 to + 0.0002406a; 2 



(14) 



For 124° < lo < 180° : 

/ a (w) = 0.025 + 0.015 sin(2.25 lo - 369.0) . 

where lo is the scattering angle in degrees. The total inte- 
grated scattering cross section iVaffa is given by eq. ( |l2"| ) for 
a given K. In B band <r a is 1.216 times the <r a in V band. 



4.3 Upward emission function 

The normalized emission function of each area gives the rel- 
ative upward flux per unit solid angle in every direction. It is 
the sum of the direct emission from fixtures and the reflected 
emission from lighted surfaces, normalized to its integral and 
is not known. The high number of sources contributing to 
the sky brightness of a site with casual distribution and ori- 
entation, smooth the shape of the average normalized emis- 
sion function which can be considered in first approximation 
axisymmetric. It is possible to measure the average upward 
emission of a chosen area at a number of different elevation 
angles when a large number of satellite measurements from 
very different orbits is available. It will be possible to obtain 
it directly by integrating upward emission from all lighting 
fixtures and all lighted surfaces on the basis of lighting en- 
gineering data and models as soon they will be available. 

In this paper we assumed that all land areas have the 
same average normalized emission function. This is equiva- 
lent to assuming that lighting habits are similar on average 
in each land area and that differences from the average are 
casually distributed in the territory. The average normalized 
emission function can be constrained from radiance mea- 
surements of cities at various distance from satellite nadir. It 
can also be obtained from comparison of Earth-based obser- 
vations and models predictions (Cinzano, in prep.). We chose 
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to assume this function and check its consistency with satel- 
lite measurements, rather than directly determine it from 
satellite measurements because at very low elevation angles 
the spread is too much large to constrain adequately the 
function shape. 

We adopted for the average normalized emission func- 
t ion t he normalized city emission function from Garstang 
(h986|) : 



= J- [2oi cos (i>) + 0.554a 2 V' 4 l 



(15) 



where ai and 02 are shape parameters. Here we assumed 
ai = 0.46, a,2 = 0.54 for the typical average function, corre- 
sponding to Garstang parameters G = 0.15 and F = 0.15. 
This function was tested by Garstang with many compar- 
isons between model predictions and measurements. This 
author assumed this function be produced by the sum of di- 
rect emission from fixtures at high zenith distances and Lam- 
bertian emission from lighted horizontal surfaces at higher 
zenith angles. Nevertheless, upward flux can be emitted at 
all zenith angles both from fixtures and vertical or horizontal 
surfaces, so we preferred to consider Garstang's function like 
a parametric representation with a\ and <i2 as shape param- 
eters without any meaning of fraction of direct and reflected 
light. We tested also the normalized city emission function 
of Cinzano (2000b, 2000c) which assumes a slightly higher 
emission at intermediate angles in respect to the function 



0.554a 2 V 



(16) 



assuming a\ = 0.79 and 0,2 = 0.21. Co mpariso n between 
these functions are presented by Cinzano (2000b). 

We checked these functions by studying the relation be- 
tween the upward flux per unit solid angle per inhabitant of 
a large number of cities and their distance from the satellite 
nadir in a single orbit satellite image taken the 13th January 
1997 h20:27 from satellite F12 which was chosen for its large 
cloudfree area. Taking into account the average orbital dis- 
tance Rs and the Earth curvature radius Rt, it is possible, 
with some geometry, to relate the distance D from satellite 
nadir with the emission angle ip: 



D . R T sm(D/R T ) 

ip = — — arcsin 

Rt s 

with: 



(17) 



8= JR 2 T + {R T + Rs) 2 



2R t (Rt +Rs) cos(-^-). (18) 
Ht 



The flux per unit solid angle per inhabitant in relative units 
is obtained dividing the measurements for the correspondant 
extinction coefficient £3(1/)) c omputed for curved-Earth from 
eq. (20) of Garstang (|1989a|): 



£3 = exp [-N m a m (Al + 11.778KA2)] 
Al = isecV> (l - e - cscos ^ + ^^^Bl 



(19) 



Bl = (c V cos 2 ip + 2cs cos ip + 2) e' cs cos * - 2 
A2 = 1-secip ( 1 _e- ascos * + ll t -^±B2] 
B2 = (a 2 s 2 cos 2 ip + 2as cos ip + 2) e^ 8008 * - 2. 

Fig. [j] shows the relative flux per unit solid angle per 
inhabitant averaged over ranges of 100 km. 
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Figure 4. Relative flux per unit solid angle per inhabitant aver- 
aged over ranges of 100 km measured for a sample of cities versus 
their distance from satellite nadir. 



100' 90 




Figure 5. Average upward flux per unit solid angle per inhab- 
itant compared with Garstang upward emission function (solid 
line) and Cinzano upward emission function (dashed line) for the 
adopted shape parameters. The components of the Garstang func- 
tion are also shown (dotted and short-dashed lines). 



In Fig. ^ we show the I(ip) obtained with eq. (Tl^) and 
( p"9| ) compared with the Garstang (1986) function (solid line) 
and with the Cinzano (20001) function (dashed line) with 
the assumed shape parameters. The fits are good with both. 
Errorbars are not necessarily related to fluctuations in func- 
tion shape but rather with fluctuations in the total flux per 
inhabitant. Given that the extinction effects seems to bal- 
ance the changes in measured flux at angles off from nadir 
(see Fig.^) , we did not need to correct the input images from 
single orbits for these off-nadir effects when computing the 
composite image. 

Snow reflects approximately 60 per cent of downlight 
changing the shape of the upward emission function. Lighted 
roads usually are cleaned in a few days so reflection of street 
lighting by snow on road surfaces is unlikely to be impor- 
tant, but reflection of the artificial sky light by snow on the 
rest of the land could enhance noticeably the upward flux in 
more polluted areas. Assuming roughly that 10 per cent of 
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Figure 6. Geometrical relationships. 



the upward flux be scattered downward by atmospheric par- 
ticles and molecules and that 60 per cent of it be reflected 
upward again by the snow covered terrain, the increase of 
artificial sky brightness by snow reflexion of sky light could 
reach 6 per cent. The upward flux due to snow reflection 
in some zones of Europe (our images are taken in winter) 
is likely detected by the OLS-PMT, but the different shape 
of the upward emission function could produce small errors. 
We plan to use specific satellite surveys to detect snowed ar- 
eas where the upward emission function must be corrected. 
Even the offshore lights in the North Sea, where oil and gas 
production sites are active, could have a different upward 
emission function. 



4.4 Geometric relations 

In this paper we are more interested in understanding and 
comparing light pollution distributions in the European ter- 
ritory rather than in predicting the effective sky brightness 
for observational purposes, so we computed the artificial sky 
brightness at sea level, in order to avoid introduction of alti- 
tude effects in our maps. We plan to take account of altitudes 
in a forecoming paper devoted to mapping the naked-eye 
star visibility which requires the computation of star-light 
extinction and natural sky brightness for the altitude of each 
land area. 

With the h ypothe sis of sea level, geometrical relations 
from Garstang (1989a) between quantities shown in Fig. I J 



taking into account Earth curvature, became simpler. They 
are listed in the appendix [a|. In eq. ( ft^ ) now is H = 0. 
Equations (0), (vm have been integrated by Garstang (1989a: 
eq. (18), (m, (22) and eq. (20), (21)). They are listed in 
appendix |A| too. 



Correction for double scattering of Garstang (1989e 
became: 



D s = 1+ iV m , cr m (ll.llisr/2 + |-) , 



(20) 



where /i and /2 are given in eq. (A3). Integration of eq. 
(^) must be done only where the Earth curvature does not 
shield the line of sight from the source. We need to start 
integration from uq as given by (from eq. (12) and (13) of 
Garstang 1989a): 



uo 



2R T sm 2 (D/2R 7 



sin z cos j3 sin(D / Rt) + 



cos z cos 



■s(D/R T ) 



(21) 



We neglected the presence of mountains which might 
shield the light emitted from the sources from a fraction 
of the atmospheric particles along the line-of-sight of the 
observer. Assuming flat Earth, mountains between source 
and observatory shield the light emitted from the source 
with an angle less than 8 = arctg^ where H is the height 
of the mountain and p its distance from the source. The 
ratio between the luminance in the shielded and not shielded 
cases is given, in first approximation, by the ratio between 
the number of particles illuminated in the two cases: 



bs_ 

fans 



CT a /a(u;) [n q N a (h)dh + (T m fm(u) fng N m (h)dh 
p p 

o-a/a(w) J °° N*(h)dh + a m / m H J °° Nm(h)dh 



,(22) 



where q is the distance of the site from the source, N m (h) 
and N a (h) are respectively the vertical number densities of 
molecules and aerosols at the altitude h, a m and cr a are their 
scattering sections and roughly u ~ -n — 8. This expression 
shows that, given the vertical extent of the atmosphere in 
respect to the highness of the mountains, the shielding is 
not negligible only when the source is very near the moun- 
tain and both are quite far from the site (Garstang 1989a; 
see also Cinzano 2000c): ^ >> where ft a is the vertical 
scale height of the atmospheric aerosols. This condition in 
the considered territories usually applies to poorly lighted 
areas only, which produce little pollution. Earth curvature 
emphasizes this behaviour. 



4.5 Relation with atmospheric conditions 

The adopted modelling technique allows us to assess the 
atmospheric conditions for which a map is computed giving 
observable quantities like the vertical extinction at sea level 
in magnitudes (Garstang 1989a): 



Am = 1.0857iVn 



,O0n 



where for V band A = 550 and for B band A = 440. Ne- 
glecting Earth curvature, the horizontal visibility defined as 
the distance at which a black object would show a bright- 
ness of 0.98 of the brightness of background horizon due to 
scattered light is (Garstang 1986): 



Ax = 



3.91 

NmfiCJ-n 



+ 11.778* (^°) 



(24) 



Effects of curvature have been discussed by the cited au- 
thor. Other relations exist with measurables like the Linke 
turbidity factor for total solar radiation received on Earth 
(Garstang 1988). The optical thickness r is (from eq. (22) 
of Garstang 1986): 



r = Am/1.0857. 



(25) 
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With K=l we obtain r = 0.3 so double scattering approxi- 
mation is adequate in our map computations. In order that 
r < 0.5 we need K <2.2. 



4.6 Spectral emission 

DMSP satellites do not carry a spectrograph able to obtain 
spectra of the upward light emitted by each land area. How- 
ever it is possible to recover information about the spectral 
emission from the differential effects of extinction on light of 
different wavelengths. Calling 7(A) the specific radiance of 
each land area at wavelength A, T(A) the PMT sensitivity 
curve and £4 (A, K) the vertical extinction for a given clarity 
parameter K, the radiance r(A) measured by the OLS-PMT 
is: 



r(K) 



T(X)I(X)£4\,K) dX. 



The vertical extinction can be obtained from eq. 
£4 (A, K) = exp (piA -4 ) • exp (p 2 A'A _1 ) 

Pi = — A'm.OCmAo/c 

P2 = —11.778 N m>0 a m Xo/a, 



(26) 



(27) 



where the molecular scattering cross section <r m at Ao = 550 
nm, the molecular density at sea level N^fi and the in- 
verse scale heights of molecules and aerosols c and a are 
given in section ^l]. Equation ( p6j ) could be inverted with 
a Richardson-Lucy algorithm in order to recover I(X) from 
r(K), a function which can be obtained comparing many ra- 
diance calibrated images obtained in different atmospheric 
conditions. The average clarity parameter K for each im- 
age can be obtained comparing the variations of measured 
radiance with increasing distance from satellite nadir. 

Fig. (?) shows in the lower panel the r(K) curves pro- 
duced by the example spectra I(X) in the upper panel. 
Curves have been scaled so that r(K—l) = 1. Spectra with a 
similar effective wavelength after convolution with the PMT 
sensitivity curve, like e.g. a constant (dashed curve) , a black- 
body at 4000° K (dotted curve) and a narrow gaussian cen- 
tered at 650 nm (long dashed curve), give almost the same 
r(K) curve. However spectra with lower or higher effective 
wavelength, like a gaussian centered at 440 nm (dot-long 
dashed curve), 550 nm (solid curve) or at 800 nm (dot- 
dashed curve), give different r(K) curves. 

We obtained the maps for B and V photometric astro- 
nomical bands. The relative B-V color index of each land 
area could be inferred from the r(K) curves but requires 
some assumptions on the shape of the spectra. In order to 
be simple in this paper we assumed it to be constant ev- 
erywhere and we only took into account the different prop- 
agation of the light in the atmosphere in the two bands. We 
plan to study differences in city spectra in future works. 



4.7 Calibration 

We calibrated the maps on the basis of both (i) accurate 
measurements of sky brightness together with extinction 
from the Earth-surface and (ii) analysis of before-fly radi- 
ance calibration of OLS-PMT. 
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Figure 7. The r(K) curves (lower panel) for some example spec- 
tra (upper panel). Spectra with a similar effective wavelength 
after convolution for the PMT sensitivity curve, like e.g. a con- 
stant (dashed curve), a blackbody at 4000° K (dotted curve) and a 
narrow gaussian centered at 650 nm (long dashed curve), give al- 
most the same r(K) curve. However spectra with lower or higher 
effective wavelength, like a gaussian centered at 550 nm (solid 
curve) or at 800 nm (dot-dashed curve), give r(K) curves with 
respectively more or less pendency. 



4-7.1 Calibration with Earth-based measurements 

A detailed calibration requires sky brightness measurements 
at a large number of sites on sea level, taken in nights with 
the same vertical extinction and horizontal visibility of the 
map under calibration, averaged over many nights in order 
to smooth atmospheric fluctuations. Observations need to 
be under the atmosphere, i.e. as actually observed from the 
ground without any extinction correction applied. To ob- 
tain the artificial sky brightness it would be necessary to 



© 1994 RAS, MNRAS 000, 



10 P. Cinzano, F. Falchi, C. D. Elvidge and K. E. Baugh 



measure the natural sky brightness in some sites where the 
maps suggest that the artificial one is negligible, and sub- 
tract the mean from all measurements. Moreover given the 
fast growth rate of artificial sky brightness, which e.g. in 
Italy reaches 10 per cent per year (Cinzano 2000d), mea- 
surements have to be taken in the same period and the cal- 
ibration will refer to that time. 

Measurements of sky brightness in Europe in V and 
B bands in the period 1996-1999 are scarce, so we cali- 
brated our maps with all available measurements in the 
chosen bands taken in 1998 and in 1999 in clean or pho- 
tometric nights even if extinction was not available or not 
exactly the required one (Catanzaro & Catalano 2000; Cin- 
zano 2000d; Delia Prugna 1999; Favero et al. 2000; Piersi- 
moni, Di Paolantonio & Brocato 2000; Poretti & Scardia 
2000; Zitelli 2000). Some data has been taken by one of us 
for this purpose with a small portable telescope and a CCD 
device (Falchi 1999). Most of the sites are at sea level but 
we included also a few sites at altitudes under 1300 m o. 
s. 1.. In lack of measurements of natural sky brightness in 
Europe at sea level, we assumed it at minimum solar activ- 
ity B = 22.7 mag arcsec -1 in B band, and V = 21.6 mag 
arcsec -1 in V band, estimating an incertitude of at least 
±0.1 mag arcsec -1 . Natural sky brightness increases when 
solar activity increases (Walker 1988) and the solar activity 
in 1998 was close to minimum but not at the minimum, so 
it could be underestimated and consequently the artificial 
brightness in darker sites considerably overestimated. The 
sky brightness has been transformed into photon radiance 
with formulae of Garstang (1989a: eq. (28) and eq. (39)). 

A least square fit of a straight line y — a + x over 
the logarithmic measured radiances versus the logarithmic 
predicted radiances gives the logarithmic calibration coeffi- 
cients a B = -0.63 ±0.04 and a v = 0.00 ±0.04. We assumed 
unavailable the incertitudes of measurements as given by 
atmospheric conditions and emission function fluctuations. 
The uncertainty of the calibration coefficients produces an 
incertitude of about 10 per cent in the calibrated predicted 
radiances. However single data points show differences even 
of about 60 per cent, so that we consider safer to assume 
this last value as estimate of the uncertainty of our cali- 
brated predictions for a given site. More precise calibrations 
will be possible when a large number of measurements of 
sky brightness at sea level together with extinction become 
available. A large CCD measurement campaign is being set 
up. 

In Figs ^ and ^ we compared our calibrated predictions 
with the available measurements of artificial sky brightness 
respectively in V and in B bands. Photon radiance in V band 
is expressed in units of 0.3419 10 10 ph s -1 m -2 sr -1 , cor- 
responding approximately to a luminance of one /^cd m -2 , 
and in B band in units of 10 10 ph s -1 m -2 sr -1 . Errorbars 
indicate measurement errors which are much smaller than 
the effects of fluctuations in atmospheric conditions. 

We checked the effects on predictions of Figs ^ and ^ 
of changes of (i) atmospheric conditions and (ii) the shape 
of An increase of the aerosol content parameter K in- 

creases exponentially the extinction and linearly the scat- 
tering along the line of sight, so that, excluding a small 
area very near the source where the extinction is negligible, 
the predicted artificial sky brightness decreases everywhere: 
a well known result (Garstang 1986, 1989a, 2000; Cinzano 




2 3 4 
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Figure 8. Measurements of artificial sky brightness versus cali- 
brated map predictions in V band. 
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Figure 9. Measurements of artificial sky brightness versus cali- 
brated map predictions in B band. 

2000c). The decrease is larger for areas farther from sources. 
These usually are darker, so in the log-log diagram, increas- 
ing haze darker sites move toward low predicted brightness 
more than less dark sites. A change of aerosol scale length 
has a similar effect. A change of the shape of the upward 
emission function I(ip) has a different effect. Decreasing the 
relative emission at low elevations (large ip) the sky bright- 
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ness is decreased approximately proportionally in almost all 
sites, excepts very near the source (see also Cinzano & Diaz 
Castro 2000), so in a log- log diagram all sites move toward 
low predicted brightness of nearly the same quantity. This 
opens a way to detect cities with more light-wasting installa- 
tions. In fact if a series of measurements of sky brightness at 
increasing distances from a city lies on a straight line parallel 
to the calibration line but displaced toward lower predicted 
values, likely the relative emission of the city at low eleva- 
tion is higher than assumed in making the map. This is a 
typical result of poorly shielded or too much inclined street 
light. 

Shifts in measurements obtained with different instru- 
mental setups could arise because, as will be shown in Fig. 
[ic| , the average emission spectrum has typically its maxi- 
mum at a side of the V band sensitivity curve, so that the 
resulting measurements are quite sensitive to little differ- 
ences between the instrumental curve and the standard V 
band curve. 



4-7.2 Calibration with pre-fly irradiance measurements 

Map calibration based on pre-fly irradiance calibration of 
OLS PMT requires the knowledge, for each land area (i, j), 
of (i) the average vertical extinction Am during satellite 
observations and (ii) the relation between the radiance in 
the chosen photometrical band and the radiance measured 
in the PMT spectral sensitivity range, which depends on the 
emission spectra. Both of them are unknown. 

If r is the energetic radiance in [l0 10 Wcm _2 sr _1 ] , the 
upward energy flux in [W] in the PMT photometric band is: 



(28) 



where I(ip) is the upward emission function and A is the 
surface area in km 2 of the land area (i, j): 



. ( 2-kAx „ V 

An,i) = Rt) cos(Z) . 

y,JJ 1360-60-60 J w 



(29) 



with I latitude of the land area, Rt average Earth radius in 
km and Ax pixel size in arcsec. 

The photon radiance in the V photometric band corre- 
sponding to the energetic radiance measured by the PMT 
is: 



fr(i,j) 



f oc 
Jo 



f 06 



r(A)S(A)* dA 



T m/w (X)I(X) ^dX T PMr (X)S(X)dX 



sr Tvwi(x) T^ dx 

J™ TiMr(A)/(X) dA 



(30) 



where TV and Tfmt are the sensitivity curves respectively of 
V band and PMT detector, 7(A) is the energy spectrum of 
the emission from the chosen land area, 5(A) is the energy 
spectrum of the pre-fly calibration source, h is the Plank 
constant and c is the velocity of light. The second ratio is 
the effective wavelength < A > of the combination of the 
sensitivity curves of the PMT and the calibration source, 
divided by he. 

In order to check the consistency of our Earth-based V- 
band calibration with pre-fly radiance calibration of PMT 
images, we obtained a tentative map calibration assuming 
for all land areas an average vertical extinction Am = 0.3 
mag V at sea level and constructing a synthetic spectrum for 
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Figure 10. The composite spectrum of upward emission (solid 
line) used for checking the consistency of Earth-based map cali- 
bration with PMT radiance calibration. We roughly assumed 50 
per cent of the total power produced by HPS lamps and 50 per 
cent by Hg vapour lamps. Dashed line shows the V band sensi- 
tivity curve. 



a typical night-time lighting. We very roughly assumed that 
50 per cent of the total power emitted by each land area be 
produced by HPS lamps (SON standard) and 50 per cent by 
Hg vapour lamps (HQI). The composite spectrum is visible 
in Fig. [l(j (solid line) together with the V band (dashed line). 

a 1.48 x 10 18 ph 
er unit power at 
>13) and eq. ©, 



The result of integration of eq. (^) is t=- 
s W _1 , much less than the photon flux 
550 nm 



2.79 x 10 18 ph s" 1 W . Eq 
taking in account the internal constants of the program and 
omitting the cos (I) already included in it, give the V-band 
logarithmic calibration coefficient av = —0.01. In spite of 
the uncertainties both in the extinction and in the average 
emission spectrum, this calibration agree very well with the 
Earth-based calibration. 



5 RESULTS 

Figs |l2|, [llj and [l4| show the artificial sky brightness 
in Europe at sea level in V and B bands. Colours corre- 
spond to ratios between the artificial sky brightness and the 
natural sky brightness of: <0.11 (black), 0.11-0.33 (blue), 
0.33-1 (green), 1-3 (yellow), 3-9 (orange), >9 (red). Original 
maps are 4800x4800 pixel images saved in 16- bit standard 
fits format with fitsio Fortran-77 routines developed by 
HEASARC at the NASA/GSFC. Images have been anal- 
ysed with FTOOLS 4.2 analysis package by HEASARC and 
with quantum IMAGE 3.6 by Quantum Image Systems. 
Maps have been computed for clean atmosphere with aerosol 
clarity K = 1, corresponding to a vertical extinction of 
Am = 0.33 mag in V band, Am = 0.56 mag in B band, 
horizontal visibility Ax = 26 km, optical depth r = 0.3. We 
limited our computations to zenith sky brightness even 
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FIGURE 



Figure 11. Artificial sky brightness at sea level in Europe in V band for aerosol content parameter K = 1. 
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FIGURE 



Figure 12. Artificial sky brightness at sea level in Europe in V band for aerosol content parameter K = 1. 
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FIGURE 



Figure 13. Artificial sky brightness at sea level in Europe in B band for aerosol content parameter K = 1. 
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FIGURE 



Figure 14. Artificial sky brightness at sea level in Europe in B band for aerosol content parameter K = 1. 
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though our method allows determinations of brightness in 
other directions. This would be useful to predict visibility in 
large territories of particular astronomical phenomena, like 
e.g. comets. A complete mapping of the art ificial brightness 
of the sky of a site, like Cinzano (2000c), using satellite 
data instead of population data is possible (Cinzano 2000 , 
in prep.). Falchi & Cinzano (200C) and Cinzano et al. (200C) 



obtained in 1998 the first maps of the artificial sky bright- 
ness from satellite data using a DMSP single-orbit image 
and replacing / in eq. ([[]) with the Treanor Law, a semi- 
empirical law which assumes a very simplified model with 
homogeneous atmosphere, vertical heights small in relation 
to the horizontal distances, scattering limited to a cone of 
small angle co-axial with the direct beam and flat Earth. Dif- 
ferences with our maps mainly arise where Earth curvature 
play a role in limiting the propagation of light pollution. Our 
study constitute the natural improvement of their seminal 
work. 

Recommendation 1 of the IAU Commission 50 (Smith 
1979) states that the increase in sky brightness at 45° ele- 
vation due to artificial light scattered from clear sky should 
not exceed 10 per cent of the lowest natural level in any part 
of the spectrum between wavelengths 3000A and 10000A. 
So this is the level over which the sky must be considered 
"polluted". The maps shows that only a few areas in Eu- 
rope are under the limit of 10 per cent at zenith and some 
of them could still result in being quite polluted at higher 
zenith distances. 



6 CONCLUSIONS 

We presented a method to map the artificial sky brightness 
in large territories in astronomical photometric bands with 
a resolution of approximately 1 km. We computed the maps 
with detailed models for the propagation in the atmosphere 
of the upward light flux measured with DMSP satellites Op- 
erational Linescan System. The use of this modelling tech- 
nique allows us to (i) assess the atmospherical conditions for 
which the maps are computed giving observable quantities 
and (ii) take into account Earth curvature. This cannot be 
done properly when using semiempirical propagation laws. 
The use of satellite data constitutes an improvement over 
the use of population data to estimate upward flux because 
(i) it allows spatially explicit detail of the geographic distri- 
bution of emissions and (ii) some polluting sources, like e.g. 
industrial areas, ports and airports, are not well represented 
in population data. 

We presented, as an application of the described 
method, the maps of artificial sky brightness in Europe at 
sea level in V and B bands. We are extending the maps to 
the rest of the World in a fore coming World Atlas of Arti- 
ficial Sky Brightness and preparing predictions for the state 
of the night sky in future years. 
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APPENDIX A: GEOMETRICAL RELATIONS 

With the hypote sis of sea level, geometrical relations from 
Garstang (1989a) between quantities in Fig. ^became: 



s = yj(u - I) 2 + 4ul sm 2 (D/2R T ) 
with 1 = a/4R2,sui 2 (D/2Rt) 



h = Rj 



ii 2 -\-2uRt cos(z) 



oj = 9 + tp 
with 8 — arccos( 



(Al) 



Y^if si n 2(D/2R T ) 

<f> = arccos( ' + 2 ls ~" ) 
ip — arccos( — ) 



qi = i?r(sin(Z) / Rt) sin(z) cos(/3) + cos(D / Rt) cos(z)) 

q 2 = u sm{z) cos(/3) sin(D/R T ) + <?3 

q 3 = u cos(z) cos(D/R T )-2R T sin 2 (D/2R T ) 

where Rt is the Earth curvature radius. 

Equations (JiJ), (Q) have been integrated by Garstang 
(1989a: eq. (18), (19), (22) and eq. (20), (21)). For sea level 
and 2 < 90° and ip < 90°, they are: 



: exp(-jV mi ocr m (pi + 11.778-Kp 2 )) 

■ c 

: a -1 sec z (1 — exp(— au cos z 

■ (c 2 u 2 cos 2 z + 2cu cos 2 + 2) exp(— cucos 2) — 2 



£1 = 
pi 

P2 
P3 

P4 — (a^u 2, cos 2 2 + 2aitcos2 + 2) exp(— em cos 2) — 2 



1 sec 2 (1 - exp(-cu cos 2) + 

16p4 tan 2 z 
9-n-2aR T 



(A2) 



: exp(-Ar mi0 (T m (/i + 11.778iff 2 )) 

: c" 1 sec i> (1 - exp(- CS cos V) + 16 /^'* 

= o-i sec^ (1 - exp(-a S cos^) + 



(A3) 



6 
h 
h 

fi = (c 2 s 2 cos 2 ip + 2cs cos i/j + 2) exp(— cs cos tj)) — 2 
fi = (a 2 s 2 cos 2 if) + 2as cos ip + 2) exp(— as cos ip) — 2 

For tp = 90° the reader is referred to the cited paper. 

APPENDIX B: USED DMSP ORBITS 

Orbits of the DMSP satellite F12 used in the composite 
radiance calibrated image. 



199603191827 
199603192009 
199603192151 
199603211803 
199603211945 
199603212127 
199603231738 
199603231920 
199603232102 
199701061828 
199701062010 
199701062152 
199701081804 
199701081946 
199701082128 
199701101740 
199701101922 
199701102104 
199701121715 
199701121857 
199701122039 
199701122221 
199702041740 
199702041922 
199702042104 
199702061716 



199702061858 
199702062040 
199702062222 
199702081833 
199702082015 
199702082157 
199702101809 
199702101951 
199702102133 
199603171709 
199603171851 
199603172215 
199603181839 
199603182021 
199603182203 
199603201815 
199603201957 
199603202139 
199603221750 
199603221932 
199603222114 
199701051841 
199701052023 
199701052205 
199701071816 
199701072140 



199701091752 
199701091934 
199701092116 
199701111727 
199701111909 
199701112051 
199701131703 
199701131845 
199701132027 
199701132209 
199702031752 
199702031934 
199702032116 
199702032258 
199702051728 
199702051910 
199702052052 
199702071703 
199702071845 
199702072027 
199702072209 
199702091821 
199702092003 
199702092145 
199701071958 
199701091752 
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